require(readstata13)
require(mlogit)
require(MASS)

#### prepare data ####
voter.data <- read.dta13("JESV��2�g2012�N�O�@�I����X������.dta", convert.factors  =  FALSE)

## independent variables
# respondent ID
id <- voter.data[, 1]
# gender (female = 1)
gender <- voter.data$f1 - 1
# age
age <- voter.data$f2
# education
education <- voter.data$f4
# exremity
extremity <- ifelse(is.na(voter.data$q17), 0, abs(voter.data$q17 - 6))
# partisanship
partisan.D <- ifelse(is.na(voter.data$q54), 0, ifelse(voter.data$q54 == 2, 1, 0))
partisan.R <- ifelse(is.na(voter.data$q54), 0, ifelse(voter.data$q54 == 11, 1, 0))
partisan.L <- ifelse(is.na(voter.data$q54), 0, ifelse(voter.data$q54 == 1, 1, 0))
partisan.N <- ifelse(is.na(voter.data$q54), 0, ifelse(voter.data$q54 == 15, 1, 0))

## dependent variable (positioning pattern)
issue.names <- c("fiscal reform", "collective self-defense", "welfare", 
                 "nuclear power", "decentralization", "constitution", 
                 "public pensions", "Futenma military base", "free competition")
issue.labels <- c("eco", "csd", "wel", "nuk", "loc", "con", "pen", "okn", "eql")
perception <- list(eco = data.frame(voter.data[, 63:73]),  # fiscal reform
                   csd = data.frame(voter.data[, 76:86]),  # collective self-defense
                   wel = data.frame(voter.data[, 89:99]),  # welfare
                   nuk = data.frame(voter.data[, 102:112]),  # nuclear power
                   loc = data.frame(voter.data[, 115:125]),  # decentralization
                   con = data.frame(voter.data[, 128:138]),  # constitution
                   pen = data.frame(voter.data[, 141:151]),  # public pensions
                   okn = data.frame(voter.data[, 154:164]),  # Futenma military base
                   eql = data.frame(voter.data[, 167:177]))  # free competition
# party names
for(i in 1:9){
  colnames(perception[[i]]) <- c("LDP", "DPJ", "CGP", "SDP", "JCP", "PNP", 
                                 "YP", "NRP", "NPD", "TPJ", "JRP")
}
# code positioning patterns (D = D/LR, R = R/LD, L = L/DR, and N = LDR)
per.pattern <- rep(NA, nrow(voter.data))
for(i in 1:9){
  LDP.per <- ifelse(perception[[i]]$LDP < 3, 0, 1)
  DPJ.per <- ifelse(perception[[i]]$DPJ < 3, 0, 1)
  JRP.per <- ifelse(perception[[i]]$JRP < 3, 0, 1)
  pattern <- ifelse(LDP.per == DPJ.per, 
                    ifelse(LDP.per == JRP.per, "N", "R"), 
                    ifelse(DPJ.per == JRP.per, "L", "D"))
  per.pattern <- data.frame(per.pattern, 
                            factor(pattern, levels = c("D", "R", "L", "N")))
}
per.pattern <- per.pattern[, -1]
colnames(per.pattern) <- issue.labels

#### distributions of posisitoning patterns (Table 2) ####
each.ratio <- matrix(NA, 9, 4)
for (i in 1:9) {
  each.ratio[i, ] <- table(per.pattern[, i]) / sum(table(per.pattern[, i]))
}
DK.rate <- rep(NA, 9)
for (i in 1:9) {
  DK.rate[i] <- sum(is.na(per.pattern[, i])) / nrow(per.pattern)
}
distribution.patterns <- cbind(each.ratio, DK.rate)
rownames(distribution.patterns) <- issue.names
colnames(distribution.patterns) <- c("liberal-conservative", "radical-reform-maintenance", 
                                     "ambiguous", "no distinction", "don't know rate")
round(distribution.patterns, 3) * 100

#### test the central argumant (online Appendix A) ####
## number of times positioning patterns selected
per.pattern.na.omit <- na.omit(per.pattern)
pattern.sum <- matrix(NA, nrow(per.pattern.na.omit), 4)
for (i in 1:nrow(per.pattern.na.omit)) {
  pattern.sum[i, 1] <- sum(per.pattern.na.omit[i, ] == "D")
  pattern.sum[i, 2] <- sum(per.pattern.na.omit[i, ] == "R")
  pattern.sum[i, 3] <- sum(per.pattern.na.omit[i, ] == "L")
  pattern.sum[i, 4] <- sum(per.pattern.na.omit[i, ] == "N")
}
colnames(pattern.sum) <- c("D", "R", "L", "N")

##  distribution of the number of times of selection predicted by the null hypothesis
pattern.matrix <- expand.grid(x1 = c(0, 1), x2 = c(0, 1), x3 = c(0, 1), 
                              x4 = c(0, 1), x5 = c(0, 1), x6 = c(0, 1), 
                              x7 = c(0, 1), x8 = c(0, 1), x9 = c(0, 1))
null.prob <- matrix(NA, 10, 4)
for (i in 1:4) {
  each.prob <- as.matrix(pattern.matrix) * matrix(rep(each.ratio[, i], 512), 512, 9, byrow = TRUE) + 
    (1 - as.matrix(pattern.matrix)) * (1 - matrix(rep(each.ratio[, i], 512), 512, 9, byrow = TRUE))
  pattern.prob <- apply(each.prob, 1, prod)
  for (j in 0:9) {
    null.prob[j + 1, i] <- sum(pattern.prob[rowSums(pattern.matrix) == j])
  }
}

## draw Figure A1
select.table <- matrix(NA, 10, 4)
null.variance <- actual.variance <- rep(NA, 4)
for (i in 1:4) {
  select.table[, i] <- table(factor(pattern.sum[, i], levels = as.character(c(0:9))))
  actual.variance[i] <- sprintf("%4.2f", var(pattern.sum[, i]))
  null.variance[i] <- sprintf("%4.2f", 
                              sum(null.prob[, i] * c(0:9) ^ 2) - 
                                sum(null.prob[, i] * c(0:9)) ^ 2)
}

distribution.main <- c("D/LR", "R/LD", "L/DR", "LDR")
png("FigA1.png", width = 150, height = 112, units = "mm", pointsize = 8, res = 600)
layout(matrix(1:4, 2, 2, byrow = TRUE))
for (i in 1:4) {
  plot(0:9, null.prob[, i], pch = 21, ylim = c(0, 0.5), 
       main = distribution.main[i], xlab = "Times of selection", ylab = "Probability", 
       xaxt = "n", yaxt = "n")
  points(0:9, select.table[, i] / nrow(pattern.sum), pch = 19)
  legend("topright", pch = c(21, 19), legend = c(null.variance[i], actual.variance[i]), bty = "n")
  axis(1, at = 0:9)
  axis(2, at = c(0, 0.25, 0.5))
}
dev.off()

## chi-squared tests
chisq.test(select.table[, 2], p = null.prob[, 2])  # D/LR
chisq.test(select.table[, 3], p = null.prob[, 3])  # R/LD
chisq.test(select.table[, 1], p = null.prob[, 1])  # L/DR
chisq.test(select.table[, 4], p = null.prob[, 4])  # LDR
round(qchisq(0.999, df = 9))  # threshold of chi-squared values at the 0.1% level

#### estimate the mixed logit model to test Hypotheses 1-3 ####
## prepare an mlogit.data object
# variables necessary for the analysis
pattern.dataset <- NULL
for(i in 1:9){
  store <- data.frame(per.pattern = per.pattern[, i], issue = issue.labels[i], 
                      partisan.D = partisan.D, partisan.R = partisan.R, 
                      partisan.L = partisan.L, partisan.N = partisan.N, 
                      gender, age, education, extremity, id = 1:nrow(voter.data))
  pattern.dataset <- rbind(pattern.dataset, store)
}
pattern.dataset.comp <- na.omit(pattern.dataset)

# number of missing observations
independent.vars <- cbind(gender, age, education, extremity, 
                          partisan.D = partisan.D, partisan.R = partisan.R, 
                          partisan.L = partisan.L, partisan.N = partisan.N)
independent.vars <- na.omit(independent.vars)
# proportion of missing observations (p. 683, lines 14-15)
nrow(voter.data)
nrow(voter.data) - nrow(independent.vars)
round((nrow(voter.data) - nrow(independent.vars)) / nrow(independent.vars) * 100, 1)

# descriptive statistics (Table A2)
basic.statistic <- matrix(NA, 8, 4)
for (i in 1:8) {
  basic.statistic[i, 1] <- mean(independent.vars[, i])
  basic.statistic[i, 2] <- sd(independent.vars[, i])
  basic.statistic[i, 3] <- quantile(independent.vars[, i], prob = 0.1)
  basic.statistic[i, 4] <- quantile(independent.vars[, i], prob = 0.9)
}
round(basic.statistic, 2)

# convert a data frame to an mlogit.data object
dataset <- mlogit.data(pattern.dataset.comp, id.var = "id", choice = "per.pattern", 
                       shape = "wide", varying = 3:6, sep = ".")

## estimation
set.seed(05222015)
mlogit.result <- mlogit(per.pattern ~ 0 | issue + gender + age + education + extremity | partisan, 
                        data = dataset, reflevel = "N", 
                        rpar = c("D:(intercept)" = "n", "L:(intercept)" = "n", "R:(intercept)" = "n"), 
                        correlation = TRUE, R = 100, halton = NA, panel = TRUE, print.level = 1)
# save the mlogit result object
# save(mlogit.result, file = "mlogit_result.Rdata")
# load the mlogit result object
# load("mlogit_result.Rdata")

#### estimation results of parameters (Table A3) ####
## extract estimates from the mlogit result object
# coefficients
coef <- mlogit.result$coefficients
# variance-covariance matrix
vcv <- vcov(mlogit.result)
# standard errors
se <- sqrt(diag(vcv))
# p values
pval <- summary(mlogit.result)$CoefTable[, 4]
# standard deviation of disturbance terms
sd.est <- sqrt(diag(cov.mlogit(mlogit.result)))
# correlation of disturbance terms
cor.est <- as.vector(cor.mlogit(mlogit.result))[c(2, 3, 6)]

# compute uncertainty in the covariance matrix of disturbance terms
var.sim.result <- matrix(NA, 10000, 6)
for (i in 1:10000) {
  coef.sim <- mvrnorm(1, coef, vcv)[(length(coef) - 5):length(coef)]
  upper <- matrix(c(coef.sim[1:3], 0, coef.sim[4:5], 0, 0, coef.sim[6]), 3, 3, byrow = TRUE)
  vcov.sim <- t(upper) %*% upper
  sd.sim <- sqrt(diag(vcov.sim))
  var.sim.result[i, 1:3] <- sd.sim
  var.sim.result[i, 4] <- vcov.sim[1, 2] / (sd.sim[1] * sd.sim[2])
  var.sim.result[i, 5] <- vcov.sim[1, 3] / (sd.sim[1] * sd.sim[3])
  var.sim.result[i, 6] <- vcov.sim[2, 3] / (sd.sim[2] * sd.sim[3])
}
var.se <- apply(var.sim.result, 2, sd)

# function to add asterisks
significance <- function(p.value) {
  star <- ifelse(p.value < 0.01, "**", 
                 ifelse(p.value < 0.05, "*", ""))
  star
}

## display the table
coef.table <- matrix(NA, 20, 4)
rownames(coef.table) <- c("Individual-specific variables", " Female", "1", " Age", "2", " educationcation", "3", 
                          " Extremity", "4", "Alternative-specific variables", 
                          " Partisanship", "6", "Disturbance term", " Standard deviation", "7", 
                          " Correlation of eta_1k and eta_2k", " Correlation of eta_1k and eta_3k", 
                          " Correlation of eta_2k and eta_3k", 
                          "Number of responses", "Number of respondents")
colnames(coef.table) <- c("D/LR", "R/LD", "L/DR", "LDR")
coef.table[c(seq(2, 8, by = 2), 11), 1] <- paste(sprintf("%5.3f", coef[c(seq(28, 37, by = 3), 41)]), 
                                                 significance(pval[c(seq(28, 37, by = 3), 41)]), sep = "")
coef.table[c(seq(2, 8, by = 2), 11), 2] <- paste(sprintf("%5.3f", coef[c(seq(30, 39, by = 3), 43)]), 
                                                 significance(pval[c(seq(30, 39, by = 3), 43)]), sep = "")
coef.table[c(seq(2, 8, by = 2), 11), 3] <- paste(sprintf("%5.3f", coef[c(seq(29, 38, by = 3), 42)]), 
                                                 significance(pval[c(seq(29, 38, by = 3), 42)]), sep = "")
coef.table[c(seq(3, 9, by = 2), 12), 1] <- paste("(", sprintf("%5.3f", se[c(seq(28, 37, by = 3), 41)]), ")", sep = "")
coef.table[c(seq(3, 9, by = 2), 12), 2] <- paste("(", sprintf("%5.3f", se[c(seq(30, 39, by = 3), 43)]), ")", sep = "")
coef.table[c(seq(3, 9, by = 2), 12), 3] <- paste("(", sprintf("%5.3f", se[c(seq(29, 38, by = 3), 42)]), ")", sep = "")
coef.table[11, 4] <- paste(sprintf("%5.3f", coef[40]), significance(pval[40]), sep = "")
coef.table[12, 4] <- paste("(", sprintf("%5.3f", se[40]), ")", sep = "")
coef.table[14, 1:3] <- sprintf("%5.3f", sd.est[c(1, 3, 2)])
coef.table[15, 1:3] <- paste(" (", sprintf("%5.3f", var.se[c(1, 3, 2)]), ")", sep = "")
coef.table[16:18, 2] <- paste(sprintf("%5.3f", cor.est[c(2, 1, 3)]), 
                              " (", sprintf("%5.3f", var.se[c(5, 4, 6)]), ")", sep = "")
coef.table[19, 2] <- nrow(pattern.dataset.comp)
coef.table[20, 2] <- nlevels(as.factor(pattern.dataset.comp$id))
print(coef.table, quote = FALSE)

#### simulation of the probability of the positioning patterns (Figures 3 and A3) ####
## prepare functions
# function to simulate probability for individual-specific variables
prob.sim.ind <- function(issue, ind.var0, ind.var1, alt.var, sim, seed) {
  set.seed(seed)
  coef.sample <- mvrnorm(sim, coef, vcv)
  if (issue == 1) {
    intercept <- coef.sample[, 1:3]
  } else {
    intercept <- coef.sample[, 1:3] + coef.sample[, (3 * issue - 2):(3 * issue)]
  }
  par1 <- cbind(intercept[, 1], coef.sample[, seq(28, 37, by = 3)])
  par2 <- cbind(intercept[, 3], coef.sample[, seq(30, 39, by = 3)])
  par3 <- cbind(intercept[, 2], coef.sample[, seq(29, 38, by = 3)])
  xb1.0 <- exp(ind.var0 %*% t(par1) + alt.var[, 2] %*% t(coef.sample[, 41])) 
  xb2.0 <- exp(ind.var0 %*% t(par2) + alt.var[, 3] %*% t(coef.sample[, 43]))
  xb3.0 <- exp(ind.var0 %*% t(par3) + alt.var[, 1] %*% t(coef.sample[, 42]))
  xb4 <- exp(alt.var[, 4] %*% t(coef.sample[, 40]))
  xb1.1 <- exp(ind.var1 %*% t(par1) + alt.var[, 2] %*% t(coef.sample[, 41]))
  xb2.1 <- exp(ind.var1 %*% t(par2) + alt.var[, 3] %*% t(coef.sample[, 43]))
  xb3.1 <- exp(ind.var1 %*% t(par3) + alt.var[, 1] %*% t(coef.sample[, 42]))
  denom0 <- xb1.0 + xb2.0 + xb3.0 + xb4
  denom1 <- xb1.1 + xb2.1 + xb3.1 + xb4
  prob1 <- colMeans(xb1.1 / denom1 - xb1.0 / denom0)
  prob2 <- colMeans(xb2.1 / denom1 - xb2.0 / denom0)
  prob3 <- colMeans(xb3.1 / denom1 - xb3.0 / denom0)
  prob4 <- colMeans(xb4 / denom1 - xb4 / denom0)
  result <- cbind(c(quantile(prob1, probs = 0.025), quantile(prob2, probs = 0.025), 
                    quantile(prob3, probs = 0.025), quantile(prob4, probs = 0.025)), 
                  c(quantile(prob1, probs = 0.05), quantile(prob2, probs = 0.05), 
                    quantile(prob3, probs = 0.05), quantile(prob4, probs = 0.05)), 
                  c(mean(prob1), mean(prob2), mean(prob3), mean(prob4)), 
                  c(quantile(prob1, probs = 0.95), quantile(prob2, probs = 0.95), 
                    quantile(prob3, probs = 0.95), quantile(prob4, probs = 0.95)), 
                  c(quantile(prob1, probs = 0.975), quantile(prob2, probs = 0.975), 
                    quantile(prob3, probs = 0.975), quantile(prob4, probs = 0.975)))
  colnames(result) <- c("2.5%", "5%", "mean", "95%", "97.5%")
  result
}

# function to simulate probability for alternative-specific variables
prob.sim.alt <- function(issue, ind.var, alt.var0, alt.var1, sim, seed) {
  set.seed(seed)
  coef.sample <- mvrnorm(sim, coef, vcv)
  if (issue == 1) {
    intercept <- coef.sample[, 1:3]
  } else {
    intercept <- coef.sample[, 1:3] + coef.sample[, (3 * issue - 2):(3 * issue)]
  }
  par1 <- cbind(intercept[, 1], coef.sample[, seq(28, 37, by = 3)])
  par2 <- cbind(intercept[, 3], coef.sample[, seq(30, 39, by = 3)])
  par3 <- cbind(intercept[, 2], coef.sample[, seq(29, 38, by = 3)])
  xb1.0 <- exp(ind.var %*% t(par1) + alt.var0[, 2] %*% t(coef.sample[, 41])) 
  xb2.0 <- exp(ind.var %*% t(par2) + alt.var0[, 3] %*% t(coef.sample[, 43]))
  xb3.0 <- exp(ind.var %*% t(par3) + alt.var0[, 1] %*% t(coef.sample[, 42]))
  xb4.0 <- exp(alt.var0[, 4] %*% t(coef.sample[, 40]))
  xb1.1 <- exp(ind.var %*% t(par1) + alt.var1[, 2] %*% t(coef.sample[, 41]))
  xb2.1 <- exp(ind.var %*% t(par2) + alt.var1[, 3] %*% t(coef.sample[, 43]))
  xb3.1 <- exp(ind.var %*% t(par3) + alt.var1[, 1] %*% t(coef.sample[, 42]))
  xb4.1 <- exp(alt.var1[, 4] %*% t(coef.sample[, 40]))
  denom0 <- xb1.0 + xb2.0 + xb3.0 + xb4.0
  denom1 <- xb1.1 + xb2.1 + xb3.1 + xb4.1
  prob1 <- colMeans(xb1.1 / denom1 - xb1.0 / denom0)
  prob2 <- colMeans(xb2.1 / denom1 - xb2.0 / denom0)
  prob3 <- colMeans(xb3.1 / denom1 - xb3.0 / denom0)
  prob4 <- colMeans(xb4.1 / denom1 - xb4.0 / denom0)
  result <- cbind(c(quantile(prob1, probs = 0.025), quantile(prob2, probs = 0.025), 
                    quantile(prob3, probs = 0.025), quantile(prob4, probs = 0.025)), 
                  c(quantile(prob1, probs = 0.05), quantile(prob2, probs = 0.05), 
                    quantile(prob3, probs = 0.05), quantile(prob4, probs = 0.05)), 
                  c(mean(prob1), mean(prob2), mean(prob3), mean(prob4)), 
                  c(quantile(prob1, probs = 0.95), quantile(prob2, probs = 0.95), 
                    quantile(prob3, probs = 0.95), quantile(prob4, probs = 0.95)), 
                  c(quantile(prob1, probs = 0.975), quantile(prob2, probs = 0.975), 
                    quantile(prob3, probs = 0.975), quantile(prob4, probs = 0.975)))
  colnames(result) <- c("2.5%", "5%", "mean", "95%", "97.5%")
  result
}

# values used for simulations (10 and 90 percentiles)
sim.value <- matrix(NA, 4, 2)
for (i in 1:4) {
  sim.value[i, 1] <- quantile(independent.vars[, i], 0.1)
  sim.value[i, 2] <- quantile(independent.vars[, i], 0.9)
}

## conduct simulations
# collective self-defense
alt.var <- subset(pattern.dataset.comp, issue  == "csd")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "csd")[, 7:10]))
sim.result.1 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "csd")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.1[[i]] <- prob.sim.ind(2, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.1[[i + 4]] <- prob.sim.alt(2, ind.var, alt.var0, alt.var1, 10000, 11192014)
}
names(sim.result.1) <- c("gender", "age", "educationcation", "extremity", 
                         "LDP partisans", "DPJ partisans", "JRP partisans", "independent")
sim.result.1

# decentralization
alt.var <- subset(pattern.dataset.comp, issue  == "loc")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "loc")[, 7:10]))
sim.result.2 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "loc")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.2[[i]] <- prob.sim.ind(5, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.2[[i + 4]] <- prob.sim.alt(5, ind.var, alt.var0, alt.var1, 10000, 11192014)
}
names(sim.result.2) <- c("gender", "age", "educationcation", "extremity", 
                         "LDP partisans", "DPJ partisans", "JRP partisans", "independent")
sim.result.2

# fiscal reform
alt.var <- subset(pattern.dataset.comp, issue  == "eco")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "eco")[, 7:10]))
sim.result.3 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "eco")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.3[[i]] <- prob.sim.ind(1, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.3[[i + 4]] <- prob.sim.alt(1, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

# welfare
alt.var <- subset(pattern.dataset.comp, issue  == "wel")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "wel")[, 7:10]))
sim.result.4 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "wel")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.4[[i]] <- prob.sim.ind(3, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.4[[i + 4]] <- prob.sim.alt(3, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

# nuclear power
alt.var <- subset(pattern.dataset.comp, issue  == "nuk")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "nuk")[, 7:10]))
sim.result.5 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "nuk")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.5[[i]] <- prob.sim.ind(4, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.5[[i + 4]] <- prob.sim.alt(4, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

# constitution
alt.var <- subset(pattern.dataset.comp, issue  == "con")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "con")[, 7:10]))
sim.result.6 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "con")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.6[[i]] <- prob.sim.ind(6, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.6[[i + 4]] <- prob.sim.alt(6, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

# public pensions
alt.var <- subset(pattern.dataset.comp, issue  == "pen")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "pen")[, 7:10]))
sim.result.7 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "pen")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.7[[i]] <- prob.sim.ind(7, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.7[[i + 4]] <- prob.sim.alt(7, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

# Futenma military base
alt.var <- subset(pattern.dataset.comp, issue  == "okn")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "okn")[, 7:10]))
sim.result.8 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "okn")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.8[[i]] <- prob.sim.ind(8, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.8[[i + 4]] <- prob.sim.alt(8, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

# free competition
alt.var <- subset(pattern.dataset.comp, issue  == "eql")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "eql")[, 7:10]))
sim.result.9 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp, issue  == "eql")[, 7:10]))
  ind.var0[, i + 1] <- sim.value[i, 1]
  ind.var1[, i + 1] <- sim.value[i, 2]
  sim.result.9[[i]] <- prob.sim.ind(9, ind.var0, ind.var1, alt.var, 10000, 11192014)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.9[[i + 4]] <- prob.sim.alt(9, ind.var, alt.var0, alt.var1, 10000, 11192014)
}

## draw the figures
# Figure 3
colors2 <- c("#0041ff", "#35a16b", "#ff2800", "#7f878f")
tiff("Fig3.tif",  width = 6, height = 7.2, units = "in", pointsize = 8, res = 300)
layout(matrix(c(1, 2), 1, 2))
par(mar = c(4, 4, 4, 1))
dotchart(c(rev(sim.result.1[[1]][, 3]), rev(sim.result.1[[2]][, 3]), 
           rev(sim.result.1[[3]][, 3]), rev(sim.result.1[[4]][, 3]), 
           rev(sim.result.1[[5]][, 3]), rev(sim.result.1[[6]][, 3]), 
           rev(sim.result.1[[7]][, 3]), rev(sim.result.1[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 19, xlim = c(-0.15, 0.15), main = "Collective Self-Defense")
abline(v = 0, lty = 3, col = "gray")
for (i in 1:8) {
  for (j in 1:4) {
    points(sim.result.1[[i]][j, 3], (53 - 6 * i - j), pch = 19, col = colors2[j])
    segments(sim.result.1[[i]][j, 1], (53 - 6 * i - j), sim.result.1[[i]][j, 5], (53 - 6 * i - j), lwd = 1.5, col = colors2[j])
  }
}
dotchart(c(rev(sim.result.2[[1]][, 3]), rev(sim.result.2[[2]][, 3]), 
           rev(sim.result.2[[3]][, 3]), rev(sim.result.2[[4]][, 3]), 
           rev(sim.result.2[[5]][, 3]), rev(sim.result.2[[6]][, 3]), 
           rev(sim.result.2[[7]][, 3]), rev(sim.result.2[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 19, xlim = c(-0.15, 0.15), main = "Decentralization")
abline(v = 0, lty = 3, col = "gray")
for (i in 1:8) {
  for (j in 1:4) {
    points(sim.result.2[[i]][j, 3], (53 - 6 * i - j), pch = 19, col = colors2[j])
    segments(sim.result.2[[i]][j, 1], (53 - 6 * i - j), sim.result.2[[i]][j, 5], (53 - 6 * i - j), lwd = 1.5, col = colors2[j])
  }
}
dev.off()

# Figure A3
png("FigA3_1.png", width = 150, height = 180, units = "mm", pointsize = 7, res = 600)
layout(matrix(c(1, 2), 1, 2))
dotchart(c(rev(sim.result.3[[1]][, 3]), rev(sim.result.3[[2]][, 3]), 
           rev(sim.result.3[[3]][, 3]), rev(sim.result.3[[4]][, 3]), 
           rev(sim.result.3[[5]][, 3]), rev(sim.result.3[[6]][, 3]), 
           rev(sim.result.3[[7]][, 3]), rev(sim.result.3[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Fiscal Reform")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.3[[i]][j, 1], (53 - 6 * i - j), sim.result.3[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dotchart(c(rev(sim.result.4[[1]][, 3]), rev(sim.result.4[[2]][, 3]), 
           rev(sim.result.4[[3]][, 3]), rev(sim.result.4[[4]][, 3]), 
           rev(sim.result.4[[5]][, 3]), rev(sim.result.4[[6]][, 3]), 
           rev(sim.result.4[[7]][, 3]), rev(sim.result.4[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Welfare")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.4[[i]][j, 1], (53 - 6 * i - j), sim.result.4[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dev.off()

png("FigA3_2.png", width = 150, height = 180, units = "mm", pointsize = 7, res = 600)
layout(matrix(c(1, 2), 1, 2))
dotchart(c(rev(sim.result.5[[1]][, 3]), rev(sim.result.5[[2]][, 3]), 
           rev(sim.result.5[[3]][, 3]), rev(sim.result.5[[4]][, 3]), 
           rev(sim.result.5[[5]][, 3]), rev(sim.result.5[[6]][, 3]), 
           rev(sim.result.5[[7]][, 3]), rev(sim.result.5[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Nuclear Power")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.5[[i]][j, 1], (53 - 6 * i - j), sim.result.5[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dotchart(c(rev(sim.result.6[[1]][, 3]), rev(sim.result.6[[2]][, 3]), 
           rev(sim.result.6[[3]][, 3]), rev(sim.result.6[[4]][, 3]), 
           rev(sim.result.6[[5]][, 3]), rev(sim.result.6[[6]][, 3]), 
           rev(sim.result.6[[7]][, 3]), rev(sim.result.6[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Constitution")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.6[[i]][j, 1], (53 - 6 * i - j), sim.result.6[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dev.off()

png("FigA3_3.png", width = 150, height = 180, units = "mm", pointsize = 7, res = 600)
layout(matrix(c(1, 2), 1, 2))
dotchart(c(rev(sim.result.7[[1]][, 3]), rev(sim.result.7[[2]][, 3]), 
           rev(sim.result.7[[3]][, 3]), rev(sim.result.7[[4]][, 3]), 
           rev(sim.result.7[[5]][, 3]), rev(sim.result.7[[6]][, 3]), 
           rev(sim.result.7[[7]][, 3]), rev(sim.result.7[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Public Pensions")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.7[[i]][j, 1], (53 - 6 * i - j), sim.result.7[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dotchart(c(rev(sim.result.8[[1]][, 3]), rev(sim.result.8[[2]][, 3]), 
           rev(sim.result.8[[3]][, 3]), rev(sim.result.8[[4]][, 3]), 
           rev(sim.result.8[[5]][, 3]), rev(sim.result.8[[6]][, 3]), 
           rev(sim.result.8[[7]][, 3]), rev(sim.result.8[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Futenma Military Base")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.8[[i]][j, 1], (53 - 6 * i - j), sim.result.8[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dev.off()

png("FigA3_4.png", width = 75, height = 180, units = "mm", pointsize = 7, res = 600)
dotchart(c(rev(sim.result.9[[1]][, 3]), rev(sim.result.9[[2]][, 3]), 
           rev(sim.result.9[[3]][, 3]), rev(sim.result.9[[4]][, 3]), 
           rev(sim.result.9[[5]][, 3]), rev(sim.result.9[[6]][, 3]), 
           rev(sim.result.9[[7]][, 3]), rev(sim.result.9[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Free Competition")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.9[[i]][j, 1], (53 - 6 * i - j), sim.result.9[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dev.off()

#### robustness check (online Appendix J) ####
## addutional independent variables
# employed
employed <- ifelse(voter.data$f5 == 7, NA, ifelse(voter.data$f5 == 1, 1, 0))
# independent businesses
indep.business <- ifelse(voter.data$f5 == 7, NA, ifelse(voter.data$f5 == 2 | voter.data$f5 == 3, 1, 0))
# student
student <- ifelse(voter.data$f5 == 7, NA, ifelse(voter.data$f5 == 4, 1, 0))
# homemaker
homemaker <- ifelse(voter.data$f5 == 7, NA, ifelse(voter.data$f5 == 5, 1, 0))
# own house
own.house <- ifelse(voter.data$f7 < 3, 1, 0)
# income
income <- voter.data$f8

## prepare an mlogit.data object
# variables necessary for the analysis
pattern.dataset.2 <- NULL
for(i in 1:9){
  store <- data.frame(per.pattern = per.pattern[, i], issue = issue.labels[i], 
                      partisan.L = partisan.L, partisan.D = partisan.D, 
                      partisan.R = partisan.R, partisan.N = partisan.N, 
                      gender, age, education, extremity, 
                      employed, indep.business, student, homemaker, own.house, income, 
                      id = 1:nrow(voter.data))
  pattern.dataset.2 <- rbind(pattern.dataset.2, store)
}
pattern.dataset.comp.2 <- na.omit(pattern.dataset.2)

# number of missing observations
independent.vars.2 <- cbind(gender, age, education, extremity, 
                            employed, indep.business, student, homemaker, own.house, income, 
                            partisan.D = partisan.D, partisan.R = partisan.R, 
                            partisan.L = partisan.L, partisan.N = partisan.N)
independent.vars.2 <- na.omit(independent.vars.2)
# proportion of missing observations (online Appendix p. 20)
nrow(voter.data)
nrow(voter.data) - nrow(independent.vars.2)
round((nrow(voter.data) - nrow(independent.vars.2)) / nrow(independent.vars.2) * 100, 1)

# descriptive statistics (Table A4)
basic.statistic.2 <- matrix(NA, 14, 4)
for (i in 1:14) {
  basic.statistic.2[i, 1] <- mean(independent.vars.2[, i])
  basic.statistic.2[i, 2] <- sd(independent.vars.2[, i])
  basic.statistic.2[i, 3] <- quantile(independent.vars.2[, i], prob = 0.1)
  basic.statistic.2[i, 4] <- quantile(independent.vars.2[, i], prob = 0.9)
}
round(basic.statistic.2, 2)

# convert a data frame to an mlogit.data object
dataset.2 <- mlogit.data(pattern.dataset.comp.2, id.var = "id", choice = "per.pattern", 
                         shape = "wide", varying = 3:6, sep = ".")

## estimation
set.seed(01302016)
mlogit.result.2 <- mlogit(per.pattern ~ 0 | issue + gender + age + education + extremity + 
                            employed + indep.business + student + homemaker + own.house + income| partisan, 
                          data = dataset.2, reflevel = "N", 
                          rpar = c("D:(intercept)" = "n", "L:(intercept)" = "n", "R:(intercept)" = "n"), 
                          correlation = TRUE, R = 100, halton = NA, panel = TRUE, print.level = 1)
# save the mlogit result object
# save(mlogit.result.2, file = "mlogit_result_2.Rdata")
# load the mlogit result object
# load("mlogit_result_2.Rdata")

## estimation results of parameters (Table A5)
# coefficients
coef.2 <- mlogit.result.2$coefficients
# variance-covariance matrix
vcv.2 <- vcov(mlogit.result.2)
# standard errors
se.2 <- sqrt(diag(vcv.2))
# p values
pval.2 <- summary(mlogit.result.2)$CoefTable[, 4]
# standard deviation of disturbance terms
sd.est.2 <- sqrt(diag(cov.mlogit(mlogit.result.2)))
# correlation of disturbance terms
cor.est.2 <- as.vector(cor.mlogit(mlogit.result.2))[c(2, 3, 6)]

# compute uncertainty in the covariance matrix of disturbance terms
var.sim.result <- matrix(NA, 10000, 6)
for (i in 1:10000) {  # uncertainty in the covariance matrix of disturbance terms
  coef.sim <- mvrnorm(1, coef.2, vcv.2)[(length(coef.2) - 5):length(coef.2)]
  upper <- matrix(c(coef.sim[1:3], 0, coef.sim[4:5], 0, 0, coef.sim[6]), 3, 3, byrow = TRUE)
  vcov.sim <- t(upper) %*% upper
  sd.sim <- sqrt(diag(vcov.sim))
  var.sim.result[i, 1:3] <- sd.sim
  var.sim.result[i, 4] <- vcov.sim[1, 2] / (sd.sim[1] * sd.sim[2])
  var.sim.result[i, 5] <- vcov.sim[1, 3] / (sd.sim[1] * sd.sim[3])
  var.sim.result[i, 6] <- vcov.sim[2, 3] / (sd.sim[2] * sd.sim[3])
}
var.se.2 <- apply(var.sim.result, 2, sd)

## display the table (Table A5)
coef.table.2 <- matrix(NA, 32, 4)
rownames(coef.table.2) <- c("Individual-specific variables", " Female", "1", " Age", "2", " educationcation", "3", 
                            " Extremity", "4", " Employed", "5", " Independent businesses", "6", " Student", "7", 
                            " Homemaker", "8", " Own house", "9", " Income", "10", 
                            "Alternative-specific variables", 
                            " Partisanship", "11", "Disturbance term", " Standard deviation", "12", 
                            " Correlation of eta_1k and eta_2k", " Correlation of eta_1k and eta_3k", 
                            " Correlation of eta_2k and eta_3k", 
                            "Number of responses", "Number of respondents")
colnames(coef.table.2) <- c("D/LR", "R/LD", "L/DR", "LDR")
coef.table.2[c(seq(2, 20, by = 2), 23), 1] <- paste(sprintf("%5.3f", coef.2[c(seq(28, 55, by = 3), 59)]), 
                                                    significance(pval.2[c(seq(28, 55, by = 3), 59)]), sep = "")
coef.table.2[c(seq(2, 20, by = 2), 23), 2] <- paste(sprintf("%5.3f", coef.2[c(seq(30, 57, by = 3), 61)]), 
                                                    significance(pval.2[c(seq(30, 57, by = 3), 61)]), sep = "")
coef.table.2[c(seq(2, 20, by = 2), 23), 3] <- paste(sprintf("%5.3f", coef.2[c(seq(29, 56, by = 3), 60)]), 
                                                    significance(pval.2[c(seq(29, 56, by = 3), 60)]), sep = "")
coef.table.2[c(seq(3, 21, by = 2), 24), 1] <- paste("(", sprintf("%5.3f", se.2[c(seq(28, 55, by = 3), 59)]), 
                                                    ")", sep = "")
coef.table.2[c(seq(3, 21, by = 2), 24), 2] <- paste("(", sprintf("%5.3f", se.2[c(seq(30, 57, by = 3), 61)]), 
                                                    ")", sep = "")
coef.table.2[c(seq(3, 21, by = 2), 24), 3] <- paste("(", sprintf("%5.3f", se.2[c(seq(29, 56, by = 3), 60)]), 
                                                    ")", sep = "")
coef.table.2[23, 4] <- paste(sprintf("%5.3f", coef.2[58]), significance(pval.2[58]), sep = "")
coef.table.2[24, 4] <- paste("(", sprintf("%5.3f", se.2[58]), ")", sep = "")
coef.table.2[26, 1:3] <- sprintf("%5.3f", sd.est.2[c(1, 3, 2)])
coef.table.2[27, 1:3] <- paste(" (", sprintf("%5.3f", var.se.2[c(1, 3, 2)]), ")", sep = "")
coef.table.2[28:30, 2] <- paste(sprintf("%5.3f", cor.est.2[c(2, 1, 3)]), 
                                " (", sprintf("%5.3f", var.se.2[c(5, 4, 6)]), ")", sep = "")
coef.table.2[31, 2] <- nrow(pattern.dataset.comp.2)
coef.table.2[32, 2] <- nlevels(as.factor(pattern.dataset.comp.2$id))
print(coef.table.2, quote = FALSE)

## simulation of the probability of the positioning patterns (Figures A4)
# function to simulate probability for individual-specific variables
prob.sim.ind.2 <- function(issue, ind.var0, ind.var1, alt.var, sim, seed) {
  set.seed(seed)
  coef.sample <- mvrnorm(sim, coef.2, vcv.2)
  if (issue == 1) {
    intercept <- coef.sample[, 1:3]
  } else {
    intercept <- coef.sample[, 1:3] + coef.sample[, (3 * issue - 2):(3 * issue)]
  }
  par1 <- cbind(intercept[, 1], coef.sample[, seq(28, 55, by = 3)])
  par2 <- cbind(intercept[, 3], coef.sample[, seq(30, 57, by = 3)])
  par3 <- cbind(intercept[, 2], coef.sample[, seq(29, 56, by = 3)])
  xb1.0 <- exp(ind.var0 %*% t(par1) + alt.var[, 2] %*% t(coef.sample[, 59])) 
  xb2.0 <- exp(ind.var0 %*% t(par2) + alt.var[, 3] %*% t(coef.sample[, 61]))
  xb3.0 <- exp(ind.var0 %*% t(par3) + alt.var[, 1] %*% t(coef.sample[, 60]))
  xb4 <- exp(alt.var[, 4] %*% t(coef.sample[, 58]))
  xb1.1 <- exp(ind.var1 %*% t(par1) + alt.var[, 2] %*% t(coef.sample[, 59]))
  xb2.1 <- exp(ind.var1 %*% t(par2) + alt.var[, 3] %*% t(coef.sample[, 61]))
  xb3.1 <- exp(ind.var1 %*% t(par3) + alt.var[, 1] %*% t(coef.sample[, 60]))
  denom0 <- xb1.0 + xb2.0 + xb3.0 + xb4
  denom1 <- xb1.1 + xb2.1 + xb3.1 + xb4
  prob1 <- colMeans(xb1.1 / denom1 - xb1.0 / denom0)
  prob2 <- colMeans(xb2.1 / denom1 - xb2.0 / denom0)
  prob3 <- colMeans(xb3.1 / denom1 - xb3.0 / denom0)
  prob4 <- colMeans(xb4 / denom1 - xb4 / denom0)
  result <- cbind(c(quantile(prob1, probs = 0.025), quantile(prob2, probs = 0.025), 
                    quantile(prob3, probs = 0.025), quantile(prob4, probs = 0.025)), 
                  c(quantile(prob1, probs = 0.05), quantile(prob2, probs = 0.05), 
                    quantile(prob3, probs = 0.05), quantile(prob4, probs = 0.05)), 
                  c(mean(prob1), mean(prob2), mean(prob3), mean(prob4)), 
                  c(quantile(prob1, probs = 0.95), quantile(prob2, probs = 0.95), 
                    quantile(prob3, probs = 0.95), quantile(prob4, probs = 0.95)), 
                  c(quantile(prob1, probs = 0.975), quantile(prob2, probs = 0.975), 
                    quantile(prob3, probs = 0.975), quantile(prob4, probs = 0.975)))
  colnames(result) <- c("2.5%", "5%", "mean", "95%", "97.5%")
  result
}

# function to simulate probability for alternative-specific variables
prob.sim.alt.2 <- function(issue, ind.var, alt.var0, alt.var1, sim, seed) {
  set.seed(seed)
  coef.sample <- mvrnorm(sim, coef.2, vcv.2)
  if (issue == 1) {
    intercept <- coef.sample[, 1:3]
  } else {
    intercept <- coef.sample[, 1:3] + coef.sample[, (3 * issue - 2):(3 * issue)]
  }
  par1 <- cbind(intercept[, 1], coef.sample[, seq(28, 55, by = 3)])
  par2 <- cbind(intercept[, 3], coef.sample[, seq(30, 57, by = 3)])
  par3 <- cbind(intercept[, 2], coef.sample[, seq(29, 56, by = 3)])
  xb1.0 <- exp(ind.var %*% t(par1) + alt.var0[, 2] %*% t(coef.sample[, 59])) 
  xb2.0 <- exp(ind.var %*% t(par2) + alt.var0[, 3] %*% t(coef.sample[, 61]))
  xb3.0 <- exp(ind.var %*% t(par3) + alt.var0[, 1] %*% t(coef.sample[, 60]))
  xb4.0 <- exp(alt.var0[, 4] %*% t(coef.sample[, 58]))
  xb1.1 <- exp(ind.var %*% t(par1) + alt.var1[, 2] %*% t(coef.sample[, 59]))
  xb2.1 <- exp(ind.var %*% t(par2) + alt.var1[, 3] %*% t(coef.sample[, 61]))
  xb3.1 <- exp(ind.var %*% t(par3) + alt.var1[, 1] %*% t(coef.sample[, 60]))
  xb4.1 <- exp(alt.var1[, 4] %*% t(coef.sample[, 58]))
  denom0 <- xb1.0 + xb2.0 + xb3.0 + xb4.0
  denom1 <- xb1.1 + xb2.1 + xb3.1 + xb4.1
  prob1 <- colMeans(xb1.1 / denom1 - xb1.0 / denom0)
  prob2 <- colMeans(xb2.1 / denom1 - xb2.0 / denom0)
  prob3 <- colMeans(xb3.1 / denom1 - xb3.0 / denom0)
  prob4 <- colMeans(xb4.1 / denom1 - xb4.0 / denom0)
  result <- cbind(c(quantile(prob1, probs = 0.025), quantile(prob2, probs = 0.025), 
                    quantile(prob3, probs = 0.025), quantile(prob4, probs = 0.025)), 
                  c(quantile(prob1, probs = 0.05), quantile(prob2, probs = 0.05), 
                    quantile(prob3, probs = 0.05), quantile(prob4, probs = 0.05)), 
                  c(mean(prob1), mean(prob2), mean(prob3), mean(prob4)), 
                  c(quantile(prob1, probs = 0.95), quantile(prob2, probs = 0.95), 
                    quantile(prob3, probs = 0.95), quantile(prob4, probs = 0.95)), 
                  c(quantile(prob1, probs = 0.975), quantile(prob2, probs = 0.975), 
                    quantile(prob3, probs = 0.975), quantile(prob4, probs = 0.975)))
  colnames(result) <- c("2.5%", "5%", "mean", "95%", "97.5%")
  result
}

# values used for simulations (10 and 90 percentiles)
sim.value.2 <- matrix(NA, 4, 2)  # 10 percentile and 90 percentile
for (i in 1:4) {
  sim.value.2[i, 1] <- quantile(independent.vars.2[, i], 0.1)
  sim.value.2[i, 2] <- quantile(independent.vars.2[, i], 0.9)
}

# collective self-defense
alt.var <- subset(pattern.dataset.comp.2, issue  == "csd")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp.2, issue  == "csd")[, 7:16]))
sim.result.1.2 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp.2, issue  == "csd")[, 7:16]))
  ind.var0[, i + 1] <- sim.value.2[i, 1]
  ind.var1[, i + 1] <- sim.value.2[i, 2]
  sim.result.1.2[[i]] <- prob.sim.ind.2(2, ind.var0, ind.var1, alt.var, 10000, 01302016)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.1.2[[i + 4]] <- prob.sim.alt.2(2, ind.var, alt.var0, alt.var1, 10000, 01302016)
}
names(sim.result.1.2) <- c("gender", "age", "educationcation", "extremity", 
                           "LDP partisans", "DPJ partisans", "JRP partisans", "independent")
sim.result.1.2

# decentralization
alt.var <- subset(pattern.dataset.comp.2, issue  == "loc")[, 3:6]
alt.var0 <- matrix(0, nrow(alt.var), 4)
ind.var <- as.matrix(cbind(1, subset(pattern.dataset.comp.2, issue  == "loc")[, 7:16]))
sim.result.2.2 <- list()
for (i in 1:4) {
  ind.var0 <- ind.var1 <- as.matrix(cbind(1, subset(pattern.dataset.comp.2, issue  == "loc")[, 7:16]))
  ind.var0[, i + 1] <- sim.value.2[i, 1]
  ind.var1[, i + 1] <- sim.value.2[i, 2]
  sim.result.2.2[[i]] <- prob.sim.ind.2(5, ind.var0, ind.var1, alt.var, 10000, 01302016)
}
for (i in 1:4) {
  alt.var1 <- matrix(0, nrow(alt.var), 4)
  alt.var1[, i] <- 1
  sim.result.2.2[[i + 4]] <- prob.sim.alt.2(5, ind.var, alt.var0, alt.var1, 10000, 01302016)
}
names(sim.result.2.2) <- c("gender", "age", "educationcation", "extremity", 
                           "LDP partisans", "DPJ partisans", "JRP partisans", "independent")
sim.result.2.2

## draw the figure
png("FigA4.png", width = 150, height = 180, units = "mm", pointsize = 7, res = 600)
layout(matrix(c(1, 2), 1, 2))
dotchart(c(rev(sim.result.1.2[[1]][, 3]), rev(sim.result.1.2[[2]][, 3]), 
           rev(sim.result.1.2[[3]][, 3]), rev(sim.result.1.2[[4]][, 3]), 
           rev(sim.result.1.2[[5]][, 3]), rev(sim.result.1.2[[6]][, 3]), 
           rev(sim.result.1.2[[7]][, 3]), rev(sim.result.1.2[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Collective Self-Defense")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.1.2[[i]][j, 1], (53 - 6 * i - j), sim.result.1.2[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dotchart(c(rev(sim.result.2.2[[1]][, 3]), rev(sim.result.2.2[[2]][, 3]), 
           rev(sim.result.2.2[[3]][, 3]), rev(sim.result.2.2[[4]][, 3]), 
           rev(sim.result.2.2[[5]][, 3]), rev(sim.result.2.2[[6]][, 3]), 
           rev(sim.result.2.2[[7]][, 3]), rev(sim.result.2.2[[8]][, 3])), 
         labels = rep(c("No distinction", "Ambiguous", "Reform\u2014Mainte", "Lib\u2014Con"), 9), 
         groups = factor(c(rep("Female", 4), rep("Age", 4), 
                           rep("educationcation", 4), rep("Extremity", 4), 
                           rep("LDP Partisans", 4), rep("DPJ Partisans", 4), 
                           rep("JRP Partisans", 4), rep("Independents", 4)), 
                         levels = c("Female", "Age", "educationcation", 
                                    "Extremity", "LDP Partisans", "DPJ Partisans", 
                                    "JRP Partisans", "Independents")), 
         pch = 20, xlim = c(-0.15, 0.15), main = "Decentralization")
for (i in 1:8) {
  for (j in 1:4) {
    segments(sim.result.2.2[[i]][j, 1], (53 - 6 * i - j), sim.result.2.2[[i]][j, 5], (53 - 6 * i - j), lwd = 0.5)
  }
}
abline(v = 0, lty = 3)
dev.off()